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Abstract 



Hybrid systems, and Piecewise Deterministic Markov Processes in particular, are widely 
. used to model and numerically study systems exhibiting multiple time scales in biochemical 

reaction kinetics and related areas. In this paper an almost sure convergence analysis for nu- 
merical simulation algorithms for Piecewise Deterministic Markov Processes is presented. The 
discussed numerical methods arise through discretisina a constructive method defining these 
processes. The stochastic problem of simulating the random, path-dependent jump times of 
such processes is reformulated as a hitting time problem for a system of ordinary differential 
1 equations with random threshold. Then deterministic continuous methods (methods with 

f*"^ ' dense output) are serially employed to solve these problems numerically. We show that the 

0^ ■ almost sure asymptotic convergence rate of the stochastic algorithm is identical to the order 

of the embedded deterministic method. We illustrate our theoretical findings by numerical 
examples from mathematical neuroscience were Piecewise Deterministic Markov Processes are 
. used as biophysically accurate stochastic models of neuronal membranes, stochastic simula- 

tion; hybrid algorithm; almost sure convergence; Piecewise Deterministic Markov Process. 

1 Introduction 

■ i — i ■ 

■ In recent years the number of applications of hybrid stochastic processes to model systems in 

biology, (bio-)chemical reaction kinetics and mathematical neuroscience have increased rapidly. 
These models either stem from a direct modelling approach, e.g., models of excitable biological 
membranes [U [7J [Ml HZ], or result from a multiscale approximation to more accurate particle 
models exhibiting clearly separated time scales, e.g., in biochemical reaction systems pi f!4l [TTl [28 ] 
or in gene regulatory networks [31j . The mathematically correct treatment of these is within 
the framework of General Stochastic Hybrid Systems Additional to these novel 'bioscience' 
applications for hybrid systems there are also more 'classical' applications of high interest in fields 
such as control and queueing theory, or models in financial mathematics and ecology, cf. [5J [51 
[HI Uni [20] and references therein. In this paper we focus on Piecewise Deterministic Markov 
Processes (PDMPs), which arc an important class of hybrid systems including most of the hybrid 
models considered in the literature, see[TJ[l[7J[l[3[Il[ini[I71[Ml[lIl[3ni[3I] out of the above 
mentioned studies. PDMPs are strong Markov processes that combine continuous deterministic 
time evolution and discontinuous, instantaneous, random 'jump' events. Specifically, the dynamics 
of the two components are intrinsically intertwined. On the one hand, the continuous time- 
evolution, defined by ordinary differential equations (ODEs), depends on the outcomes of discrete 
events via randomly changing parameters and, on the other hand, the probability of the discrete 
events happening, i.e., a random, instantaneous change in a parameter, depends on the time- 
evolution - the path - of the continuous variables. 
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Due to the high complexity of hybrid models, particularly in biochemical applications, they are 
studied extensively and almost exclusively by numerical means. To this end either so called 
'pseudo-exact' algorithms, i.e., simple statistically approximate algorithms obtained by an ad- 
hoc model approximation, or statistically exact algorithms are employed. Here statistically exact 
means that in principle the simulation method produces paths that arc samples of the distribution 
of the underlying stochastic process in contrast to pseudo-exact methods which produce sample 
paths from a distribution approximating the distribution of the underlying stochastic process. 
However, even statistically exact algorithms are exact only in theory due to a numerical error which 
arises inevitably in solving most systems of differential equations in actual implementations. So 
far model accuracy was primarily of interested as the main source of error and thus considerations 
regarding the numerical error are in general neglected, cf., e.g., (Tl 1141 ITT]. Yet ultimately, even 
if a theoretically exact PDMP formulation of the model is considered or the model is highly 
accurate, numerical studies are conducted by numerical approximations to the PDMPs as an 
analytic representation of the paths is in general not available. Despite the importance and 
widespread use of numerical studies a numerical analysis, in particular, a thorough analytical 
investigation of the convergence and error behaviour of algorithms to be used, is still missing. 
The aim of this article is to provide a - to the best of our knowledge - first contribution towards 
this goal. We note that the statistically exact algorithms that were introduced in the studies [1] 
and [31] fall within the class we consider, however, in the present study a convergence analysis is 
carried out. 

In particular in the present study we are interested in the convergence of numerical approximations 
to PDMPs in a pathwisc sense which corresponds to the fact that numerical simulations are 
carried out path by path. The methods we present for approximating a PDMP incorporate as 
an integral part continuous ODE methods, also called methods with dense output. Apart from 
numerically solving the deterministic inter-jump dynamics the key problem in simulating a PDMP 
is simulating the random, path-dependent jump times. As shown in Section [3] this problem can be 
reformulated and then, combined with the numerical solution of the inter-jump dynamics, yields 
a hitting time problem with random threshold. This we solve using continuous ODE methods. 
The main feature of such a method is that it does not only provide a numerical approximation 
to the exact solution at discrete grid points but provides an approximation of the whole path 
over the whole interval. Essentially, continuous methods are approximations on a discrete grid 
with an interpolation formula for the intervals between the grid points. Hence, these methods are 
naturally suited for solving hitting time problems. The main result of this paper is that numerical 
approximations of PDMPs built on continuous ODE methods conserve the order of convergence 
of the underlying continuous method. That is, if an approximation is constructed using, e.g., a 
continuous Rungc-Kutta method of order p, then also the almost sure convergence of the stochastic 
approximation to the PDMP is of order p. 

As for the PDMPs discussed in the present study we distinguish several types. The main part of 
the paper is devoted to processes with jumps occurring only in a fixed subset of the components 
of a vector- valued PDMP which are otherwise piecewise constant. Furthermore we assume that 
the jump heights are discretely distributed. We consider this particular structure for processes as 
these arise in applications in mathematical neuroscience, which initially motivated this study, and 
as a multiscalc approximation to certain chemical reaction systems. In these models fast and slow 
reactions, modelled with reaction rate equations and instantaneous random jumps, respectively, 
affect different set of reactants, however, each possesses rates also depending on the other set of 
reactants. As jumps corresponds to an instantaneous event where one or more individual reactants 
change their state, in these models jump heights are typically integer valued. The second class of 
PDMPs comprises processes for which jumps may occur in all components and these components 
need not be necessarily piecewise constant, but jumps are still discretely distributed. These 
PDMPs include chemical reaction systems where fast and slow reactions may affect the same 
type of reactants. Finally, we also consider this last general type of processes with continuously 
distributed jump heights. 

The remainder of the paper is organised as follows. In Section [2] we present a brief definition 
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of PDMPs and their construction from sequences of independent, identically distributed (i.i.d.) 
standard uniform random variables. This provides a theoretically exact simulation algorithm which 
takes the role of an exact solution for our convergence analysis. Section|3]prcscnts the approximate 
algorithms we consider in this paper and also contains the main convergence theorem. The proof 
of the convergence theorem can be found in Section |4j We extend the convergence result to other 
classes of PDMPs in Section [5] Finally, in Section [6] we present some numerical experiments using 
examples from the neuroscience literature to illustrate the theoretical findings and draw some 
conclusions for the usage and implementations of the numerical methods. 

2 PDMPs and their construction from i.i.d. sequences of 
uniform random variables 

In this section we give a brief introduction to PDMPs, introducing the objects used for their 
definition and the construction of their paths. In particular, we present Algorithm[T]which provides 
the 'exact solution' the numerical approximations we present converge to. For a general discussion 
of PDMPs we refer to the monographs [5] and, more recently, [TB] or the PhD thesis of the present 
author [25] . 

For the main part of the paper we consider a PDMP (X(t)) t6 [n,T] = (X(t, ^))te[o,T] to consist 
of two qualitatively different components, i.e., X(t) = (Y(t),9(t)) G M. d+m . For these we assume 
that Y(t) possesses continuous paths in a set D C R d and 8(t) is right continuous and piecewisc 
constant in K C R m , where K is an at most countable set. We denote by t n , n > 1, the jump 
times of the component 6(t) and by E = D x K the state space of the process. Such a PDMP 
is uniquely defined (up to versions) by its characteristic triple (/, A, /i), where the component / 
is used to define the deterministic, continuous evolution of the trajectories in between jumps and 
the components A and /i yield the stochastic jump mechanism. In detail these objects are: The 
first is a measurable map / : E — > K d such that for each 8 G K the differential equation 

y = f(y,0) (i) 

possesses a unique solution on [0, T] with values in D for all initial conditions y(0) = yo G D. For 
notational convenience we introduce the map t i— > (j)(t,xo) which denotes the unique solution to 
the system 

'v\ _ ff(y,e) 



o 



(2) 



with respect to the initial value xo = (yo,6o) 6 E. It satisfies </>(f, xq) G E for all t G [0,T] and 
all initial values xo 6 E. The system (|2|) defines the paths of the PDMP in between jumps, see 
Algorithm [T] below. Thus we can think of such a PDMP as an ODE system with a parameter 
changing randomly in a certain way now described: The non-negative, measurable map A : E — > 
R+ is the instantaneous intensity of the jumps. It defines a probability distribution R + for all 
x G E via the survivor function^ 

S(t,x) = cx p(-^ A(0(s,x)) ds^j , (3) 

which is used to define the distribution of the inter-jump times. That is, S(t, x) states the prob- 
ability that there does not occur a jump in the next time interval of length t conditional on the 
process being currently in state x. We always assume that A is path-integrable, i.e., 

T 

\(4>(s, x)) ds < oo V x G E . 



1 A survivor function t i— > S(t) S [0, 1], t > 0, of a non-negative random variable states the probability that this 
random variable is larger than t, for discussion in relation to PDMPs see [161 Sec. 4.1]. 
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Finally, fi is a Markov kernel from E into R m defining the distribution of the jump heights in 
the component 9(t) conditional on the pre-jump values. We assume that •) is a discrete 
probability measure satisfying 6), {0}) = for all x = (y,9) G E and that /i is continuous 
in x, i.e., x H> A) is continuous for every Borel set A C R m . Clearly, the at most countable 
support K x of the measure fi(x, •) has to be such that the resulting post-jump values are in K, 
i.e., for all x = (y, 9) it holds that 6 + K x C A'. 

Algorithm[T]bclow can be used to prove the existence of a unique PDMP to given characteristics in 
a constructive way |Hj ■ More importantly, this algorithm also provides a theoretically exact method 
for simulating paths of a PDMP on a probability space that supports a sequence of i.i.d. standard 
uniform random variables. Here 'exact' denotes the fact that the distribution of the process defined 
by Algorithm [1] equals the distribution of the PDMP to the triple (/, A,/x). However, the adverb 
'theoretically' should emphasise the fact that in general the algorithm is an 'impossible' algorithm 
as neither the initial value problem (I VP) © nor the implicit equation (U) arising in Algorithm 
[T]can usually be solved exactly. Hence, the algorithm cannot be employed in practice to simulate 
trajectories. Yet, it can be used as a comparison to approximate algorithms and thus, in typical 
terms of numerical analysis, plays the role of an exact solution for a convergence analysis. Now, 
let (f2, J 7 , P) be a probability space which supports a sequence of i.i.d. standard uniform random 
variables U n , neN. 

Algorithm 1. An exact simulation algorithm for a PDMP is given by: 

Step 1. Fix the initial time to = and initial condition X(0) = xq = (yoj 9q) £ E and set a jump 
counter n = 1 . 

Step 2. Simulate a uniformly distributed random variable U2n-\ and solve 

S(r,X(t n )) = U 2n _ 1 (4) 

with respect to r to obtain the waiting time until the next jump time, i.e., r = t n — 
Then for t £ t n ] set 

X{t) = 4>{t-t n -i,X{t n -x)). 

Iftn ^ T, stop at t = T . 

Step 3. Otherwise, simulate a post-jump value 9 n for the piecewise constant component accord- 
ing to the distribution of the jump heights given by fi((f>(t n — t n -i,X(t n -i)), ■) via the 
uniformly distributed random variable [/ 2 „ and set 



X(t n ) 



f 4>{tn ~ t„-l,X(t n _i)) 
^ 9n 



Step 4- Set n = n + 1 and start again with Step 2. 

Algorithm [T] is almost surely well-defined unless one of the following two cases occur: first, if 
neither a unique nor a solution at all exists to equation (01 , which can occur for PDMPs, and 
second, if trajectories possess countably many jumps before time T with positive probability. In 
this study the first case is excluded by the conditions on A in the convergence theorem and hence 
for the remainder of the paper we can assume that a unique solution to ^ always exists. The 
second case is excluded if we assume that the characteristics (/, A, fj,) are such that the PDMP is a 
regular jump process. That is, there are almost surely only finitely many jumps in any finite time 
interval. A simple condition which guarantees the regularity of a PDMP is that the jump rate A 
is bounded, see [T^] which also contain more general conditions. 

We remark that the class of PDMPs to triples (/, A, /i) as given above contains some prominent 
classes of stochastic processes as special cases. A first special case is given by ODEs with Markovian 
switching. For these processes the jump rate is given by jumps of a Poisson process that is 
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independent of the paths followed in between jumps, i.e., A = const. Such equations can be 
written as stochastic differential equations driven by a Poisson process, i.e., stochastic differential 
equations of jump type (JSDEs) without a diffusion term, cf. [21]. In general, JSDEs without a 
diffusion part are special cases of PDMPs [TB] . Secondly, PDMPs arc generalisations of continuous- 
time Markov chains. These are piecewise constant processes and thus follow a trivial evolution 
between successive jump times, i.e., / = 0. However, in contrast to switching ODEs the jump 
intensity may be state dependent, i.e., A ^ const. Such processes are most importantly used 
in stochastic models of chemical reaction kinetics and are treated numerically by the Stochastic 
Simulation Algorithm (SSA) [TT|. We note that the numerical advantage of these special cases 
is that inter-jump times can be simulated exactly by sampling from an exponential distribution, 
that is, the implicit equation (0| can be solved exactly. Accordingly there exists a large amount 
of literature on simulation methods for these processes and a numerical analysis thereof, see, for 
instance, references in the recent studies [S] [2U] . 

Example. To conclude this brief introduction we present a concrete example for a PDMP arising 
in mathematical neuroscience which is describe in more detail in Section [6] The component Y(t) 
is one-dimensional taking values in [0, E^ a ] for initial conditions therein and the component 9(t) 
is 8-dimensional where K = {6 e {0, . . .,N} 8 : J2k=i °k = N} with N e N. Thus the phase 
space is E = [0, E-^ a ] x K and the characteristics are given as follows. The family of ODEs ([TJ is 
given by 

Cy = -g Na 08 (y - #Na) - 9tV , e K 
with constants g Na , i?Na, 3l > 0j an d the stochastic dynamics are given by the jump-rate 



A(M)) 



and the transition measure /i consisting of point probabilities of the form 

M ((y, o), {(-i, i, o, o, o, o, o, of}) = \ a ;; {y) J: . 

2.1 A property of the jump height distribution 

In this section we derive a consequence of the continuity of x <— > fi(x, ■) which is important for the 
proof of the convergence result. First, we recall that a standard method used to define random 
variables from uniform distributed random variables as is needed in Algorithm Q] is the inverse 
CDF method, see, e.g., (TUJ. In particular, solving ^ in Step 2 of Algorithm [T] is a version of 
the inverse CDF method for a distribution on the positive reals. We assume that also for the 
simulation of the jump heights a version of the inverse CDF method is employed, that is, there 
exists a function H : E x (0, 1) — > K such that for all measurable A C M. m it holds 

F[H(x,U) 6 4] =ft(x,A), 

cf., e.g., [9] Corol. 23.4]. Hence, the post jump value 9 n +i in Step 3 of Algorithm Q] is given by 
6„ + H(4>(t n+ i — t n ,X(t n )), t/ 2 „ + 2)- As K is an at most countable set there exists another at most 
countable set K' containing all supports K x of the distributions /j,(x, •) for all x G E. Thus the 
Markov kernel /x can be written as a linear combination of Dirac measures 
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where pk{x) is the point probability of the jump height k £ K' conditional on the realised pre-jump 
value x. The probabilities Pk{x) may very well be - and in general are - zero for certain pairs 



5 







U 2n (u>) 



1 



Figure 1: Schematic view of the domain of H for probabilities n(X(t n — ), •) with discrete support. 
The curve marks points of discontinuity. The function H is constant 6* for arguments (x, u) left 
of the curve and piecewise constant with different values right of it. 8* is the realised nth jump 
height and we see H(y, £/2n(w)) = 9* for y close enough to X(t n — ) as the probability that Uin 
lies on a discontinuity is zero. 

(x, k). We index the events in K' from 1 to \K'\ < oo and obtain a function H by defining 

i — 1 i 

H(x,u) = fc l if ^p k] (x) <u< ^^Pkjix) . (5) 

i=i j=i 

The definition as in (J5j refers to the standard simulation method for discrete probabilities from 
standard uniform random variables by 'binning' the interval (0,1), see, e.g., |10j . However, the 
property that x t— > fi(x, A) is continuous for all measurable sets A now implies that the point 
probabilities pk{x) depend continuously on x. Thus H is continuous in the first component at a 
point x, i.e., \\m z ^ x H{z 1 u) = H(x,u), unless u = J2]=iPki( x ) f° r some i = 1, . . . , \K'\. That 
is, unless u equals a discontinuity point of u ^ H(x, u). However, the set of discontinuity points 
is at most countable and hence the probability of a standard uniform random variable U taking 
such a value is zero. Therefore it holds that we can choose almost surely a z close enough to x, 
say \z — x\ < k, where n = U) depends on x and U , such that 

\H{x,U) - H{z,U)\ = 0, (6) 

i.e., for z close enough to x the realised jump heights coincide, cf. Fig[TJ 

Finally, we remark, that the function H is in general not unique. However, choosing a different 
H just results in a different version of the same PDMP on the same probability space. Thus 
throughout the paper we assume that one such H is fixed which is used in Step 3 of Algorithm [T] 
and also in Step 3 of the numerical approximation Algorithm [2] in Section [3] and it satisfies (j6|). 

3 The main convergence theorem 

In this section we first precisely state and discuss the convergence concept we are interested in 
before we continue with constructing the numerical methods by discretising Algorithm [TJ Sub- 
sequently, at the end of the section we state the central result of this study: the almost sure 
convergence of the numerical approximations and their asymptotic order of convergence. 




Figure 2: Schematic relation of the exact PDMP's continuous component Y(t) (black line) to its 
approximation Y(t) (grey line) in phase space and the approximation errors we need to control to 
obtain path- wise convergence of Algorithm [2] 

The set-up is as in Section[2J i.e., (fi, IF, P) is a probability space supporting a sequence of i.i.d. stan- 
dard uniform random variables U n , n > 1. Further, we fix a finite time interval [0, T]. Then for 
given characteristics (/, A, fi) as in Section [2] we use the sequence of uniform random variables to 
construct a regular PDMP (X (t, ui))te[o,T] on this probability space via Algorithm [TJ Moreover, 
we assume that numerical approximations (X(t,h,uj)) t ^Q^ are defined on the same probability 
space. Here h denotes a defining parameter of the numerical approximation, such as a discreti- 
sation step size. In the following we refer to the PDMP (X(t,uj)) te [ _ T ] as the exact PDMP and 
to its numerical approximation (X(t, h, w))*e[o,Tl as the approximate PDMP. Further, t n (u)) and 
t n (h,uj), n > 1, denote the jump times of the exact PDMP and its approximation, respectively. 
Finally, the number of the exact PDMP's jumps in the time interval [0,T] is denoted by N(T,uj), 
where N(T,uj) < oo almost surely. 

In this study we are interested in the pathwise convergence of global errors of the form 

lim max \X(t n (u),u)) - X(t n (h,w), h,oj)\ = (7) 

h— >0 n=l,...,N(T,Lu) 

and 

lim \X(T,lj) -X(T,h,cj)\ = (8) 

h—tO 

for almost all ui 6 O. That is, ([7]) and ([8]) state the almost sure convergence of the approximate 
PDMP to the exact PDMP at their jump times and at the endpoint, cf. Fig. [2] It is a standard 
approach in numerical analysis, whether deterministic or stochastic, to define the error of a numer- 
ical approximation by a measure of the difference of the exact and approximate process at discrete 
time points in the approximation interval [0,T] P21 [HI HI] • Here the jump times of the PDMP 
and its approximation are used as a discretisation of the time interval. However, in general, the 
exact and approximate jump times do differ and hence yield different time discretisations for the 
exact and approximate PDMP. This is in contrast to error definitions typically used in numerics. 
Moreover, typically, the number and locations of grid points in error definitions are the same for 
each trajectory and increase in number for decreasing step size h. It is exactly the opposite in the 
error definition ([7]) as N(T,uj) varies with lu £ Q but is fixed over variations in h. 
As the numerical approximation is continuous over time almost everywhere we are in principle 
able to use a 'typical' error definition. However, for ordinary or stochastic differential equations 
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the time points the error is evaluated at arise naturally as the grid points of the discretisation used 
in the numerical method. For PDMP approximations as considered in this study this is not the 
case. In general, the methods we consider do not possess the same grid points for any trajectory of 
one PDMP to approximate and they possibly do not even possess the same number of grid points. 
Therefore, setting particular but arbitrary points apart in using them in an error definition without 
a mathematical reasoning does not seem preferable to us. However, for a PDMP all pathwisc 
numerical methods necessarily need to approximate the jump times. Therefore, relating the values 
of the exact and approximate PDMP at their respective jump times is reasonable insofar as an 
accurate approximation is expected to behave after the nth jump just as the exact PDMP behaves 
after the nth jump. 

We have already mentioned that, in general, the jump times of the exact and approximate PDMP 
differ. Therefore it is not sufficient to consider the errors (J7]) and (J5]) alone. They may become 
arbitrarily small with the exact and approximate trajectories still substantially differing as possibly 
only the values at the jump times converge but the jump times themselves remain separated. Thus, 
we additionally require from a convergent method that also the jump times of the approximate 
PDMP converge to the jump times of the exact PDMP, that is, for almost all u G 

lim max \t n (ui) — t n (h, = . (9) 

h->0n=l,...,JV(7» 

Definition 1. We say an approximation is pathwise convergent if (JT])-© hold. In addition, an 
approximation is of order p if p is the largest integer such that the asymptotic behaviour of ©-([I]) 
is 0(hP). 

In order to construct an approximation algorithm for the numerical simulation of a PDMP's 
paths we start from the theoretically exact Algorithm [TJ We obtain an approximate algorithm 
by discrctising the problems © and ((¥]) and solving them numerically. A numerical solution of 
@ can be obtained by a standard ODE method, thus substituting Step 3 with an approximation 
is straightforward. However, the delicate part is numerically solving ((4]) for which obviously a 
numerical solution of ([2]) is needed as an integrand which is integrated until an a priori unspecified 
time r, cf. the definition of the survivor function ((3]). Hence, the two problems have to be solved in 
parallel including a mechanism for detecting the time r conditional on a realisation of a standard 
uniform random variable U . For an efficient solution we transform the equation for the next 
jump time (fj| into an equivalent problem. This can then be combined with the IVP ^ yielding 
a merged formulation of the two problems which allows for a numerical solution by continuous 
ODEs methods. The resulting numerical approximation is finally given by Algorithm [2] below. 

Taking the logarithm of equation ((4]) we obtain the equivalent equation 

w(r,x)= [ X((f>{s,x))ds = -log 17. (10) 



o 



Thus w(t, x) = — log S(t, x) is the logarithm of the survivor function S. Differentiation of w with 
respect to r yields 

w(t, x) = \{4>(t, x)) . 

This in turn yields that calculating the next jump time from i.e., Step 2 in Algorithm [I] is 
equivalent to solving the IVP 




f: ! 3=( x o' ,) )' i£[o - ) (u) 



\w(0)J 

which is integrated until the component w(t) hits the threshold — log{7 2 „_ 1 . The hitting time 

r = inf {t > : w(t) = - logf/ 2n -i} 
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is equivalent to the waiting time until the next jump, where, in general, r = oo may be possible 
unless, e.g., A is bounded away from zero. 

We need to remark on two aspects at this point. Firstly, the set-up of the hitting time problem 
([lip is as in [1] . The authors therein state how it can be solved using continuous methods and also 
present an ad-hoc implementation of an event detection for this specific problem. However, as the 
authors are primarily focused on the modelling and simulation results, they suppose that (|11|) can 
be computed "up to any desired accuracy and neglect the discretisation error" (J p. 6], This may 
be a reasonable assumption for the solution of any standard ODE without jumps in the parameters 
- as numerical methods are well studied in this case -, however, additional assumptions are needed 
and have to be considered such that an event detection is possible with arbitrary accuracy, cf. |29j . 
This aspect is not discussed by the authors. An algorithm for simulating a PDMP's path repeatedly 
solves an ODE system starting in the point the last event detection yielded. Thus we can only 
assume that the algorithm produces a path with any desired accuracy if also the event detection 
locates the hitting times with any desired accuracy. To repeat the purpose of this study, it is 
precisely this point we address and present for a large class of algorithms a thorough numerical 
analysis in terms of pathwise convergence and the conditions this presupposes. 
Secondly, we note that there exist different but equivalent set-ups for the event detection problem 
(fTTj) : one as stated in [31] and a second, which is essentially analogous to that, we briefly introduce 
now. Instead of manipulating equation Q to obtain ([10[) wc differentiate the survivor function S 
with respect to r and obtain 

S(t, x) = —S(t, x) \(4>(t, x)) . 
This yields instead of (fTTj) an IVP to solve given by 



f(y,0) \ fy(o)\ , x(t 

o , 0(0) =( A ^),*e[0,oo). (12) 
-S\(y,6)J \S(0)J V J 

This system is solved until the component S(t) hits the threshold given by the random variable 
U-zn-i- Once again the hitting time is exactly the time until the next jump. Clearly the two 
systems (|11|) and ([12j) arc equivalent in the sense that the y and 9 components of their respective 
solutions coincide as do the hitting times of the respective thresholds for a given realisation of 

Uln-\- 

However, with respect to theoretical analysis the set-up ([TT]) is the most feasible and also we 
conjecture that for actual implementations it is more efficient as (jll[) is a simpler type of ODE 
system than (fT2|) . For this reason, we consider in what follows IVP (fTTj) and denote its solution 
with respect to the initial value (xojO) by ip(t,xo). Concerning the different notations introduced 
so far wc summarise 




as we make use all of these notations whenever brevity or clarity demands. 

Next, we denote by tp(t, x, h), t > 0, a continuous approximate solution to the IVP ([lip with initial 
condition (x, 0) obtained by a continuous numerical ODE method with step size h. Consistently, 
the components of i\) arc denoted by <f> and w, respectively. An example of a continuous method 
is the continuous trapezoidal rule, which applied to (p} takes the form 

y((n+l)h) = y{nh) + h\f{y{nh),e)+h\f{y{{n+l)h),d) (13) 

with n = 0,...,N-1, h = T/N and for £ G [0, 1] 

y(nh+£h) = y(nh) + h\t{2 - f(y(nh),8) + h^ 2 f(y((n+l)h),0) . (14) 
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That is, (TlU)) is the usual trapezoidal rule for ODEs supplemented with an interpolation formula 
()14p to obtain an approximation over the intervals [nh, (n+l)h\ between the discrete grid points. 
We remark that although for the numerical methods in Section [5] or the trapezoidal rule just 
described the parameter h denotes an equidistant step size this need not necessarily be the case. 
On the one hand, h may denote the maximal step size used in case of variable step size methods, 
or, on the other hand, a given error tolerance, cf. a discussion of deterministic event detection in 
[29] . Essentially, h is a defining parameter of the method which, if convergent, converges to the 
exact solution for h — > . In favour of linguistic simplicity we keep referring to h as the step size 
in the remainder of the paper. Further, to keep the presentation and notation simple we restrict 
ourselves to continuous one-step methods. The implementations in Section [6] are all based on 
continuous Runge-Kutta methods discussed in [2]. However we expect the subsequent results, i.e., 
Theorem Q] and its Corollary [TJ to remain valid in the case of continuous multi-step methods. 
We refer to [2J[T3] f° r a general discussion of continuous methods and just briefly collect properties 
of these methods which we always assume to hold when employed to solve a standard ordinary 
differential equation of the form 

y = f(t,y) (15) 

on the interval [0, T]. In order to discuss these properties we employ in this paragraph the following 
notation: We denote by y(t,x) the exact solution of (|T5t with respect to the initial condition 
?/(0, x) = x and y(t,x,h) denotes the numerical approximation with respect to the same initial 
condition and step-size h. Then, firstly, the approximate solution obtained by a continuous ODE 
method satisfies a stability-type estimate of the form 

max \y(t,x) -y(t,z,h) I < c LT \x- z\ + err(T, h) (16) 
te[o,r] 

with an error function err(T, h) that satisfies err(T, h) —> for h —> 0. Moreover, the constant 
L > and the function err depend only on the right hand side of (|15p and especially in the 
case of the right hand side of ([T5|) being Lipschitz continuous, cf. conditions of Theorem [TJ L 
can be chosen to be the Lipschitz constant of / and also err depends in this case only on the 
Lipschitz constant. Further, err does not depend on the initial condition x £ E and without loss 
of generality we may assume monotonicity for the error, i.e., err(T\,K) > err(T2,h) if T\ > T%. 
For the trapezoidal rule (|T3| . ([T4[) such a function err is given for small enough h by 

err(T, h) = c(l + L~ 1 e LT )h 2 

with some appropriate constant c > if L is a Lipschitz constant to / in (|15[) . Particularly, the 
stability condition (|16[) implies that the method is uniformly convergent on [0, T] in the sense that 

lim max \y(t, x) — y(t, x, h)\ = (17) 

h-X)t€[0,T] 

for all initial conditions x. We say that the method is convergent of order p if p is the largest 
integer such that for h — > 

max \y(t,y)-y(t,y,h)\ = 0(h p ), (18) 
te[o,T] 

in which case err(T, h) = 0(h p ). Finally, we assume that the numerical approximation satisfies a 
Lipschitz condition with respect to t on [0,T], i.e., 

\y(t,x,h)-y(s,x,h)\<C\t-s\ Vt,s€[0,T], (19) 

where the Lipschitz constant C is uniform with respect to the initial condition x. In particular, 
these conditions are satisfied by the one-step methods developed in for the IVP (fTTj) with 
properties of the right hand side as specified in Theorem [T] 

Thus we obtain an approximation (X(t)) te [o,r] for the PDMP (X(i)) te [ 0) T] using a continuous 
ODE method by the following algorithm. Here and subsequently we usually omit to denote the 
dependence of the numerical approximation on the the step size h. 
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Algorithm 2. An algorithm simulating an approximation to a PDMP is given by: 



Step 1. Fix the initial time to = and initial condition X(0) = x = (y a , 9q) € E and set a jump 
counter n = 1 . 

Step 2. Simulate a uniformly distributed random variable V^n-i an d solve the IVP with 
initial condition (X(t n _i),0) numerically until 

? = inf {t > : w(t) = - log U] 

We take t n = t n -\ +t as the numerical approximation of the next jump time. Hence, we 
set 

X(t) = for * G (tn-lX). 

If tn > r, s£op at t = T . 

Step 3. Otherwise, simulate a post-jump value n for the piecewise constant component accord- 
ing to the distribution of the jump heights given by [i((j)(t n — X(t„_x)), •) via the 
uniformly distributed random variable Uin and set 

X(t ) — ( ^ n ~ t n -i, X(t n -i)) \ 



On 



Step 4- Set n = n + 1 and start again at Step 2. 



Just as hybrid models mix stochastic jump models with continuous deterministic models, hy- 
brid algorithms are essentially constructed by combining a simulation algorithm for the discrete 
stochastic events, i.e., the SSA, and a numerical method for solving differential equations. In the 
SSA terminology Algorithm [2] simulates the stochastic events by Gillespie's direct method, hence 
such an algorithm may also be called direct hybrid method, cf. [Tj. The above kind of algorithm, 
though arising naturally from the model problem, is also related to numerical methods developed 
for JSDEs. In the case of A = const. Algorithm [2] is exact and turns into jump-adapted methods 
for JSDEs, see, e.g., [5S] for a discussion of jump-adapted Taylor methods, in particular for the 
pure jump case. We note that derivation or analysis of numerical methods for PDMPs along the 
lines of Ito- Taylor expansions employed in the JSDE case is not possible as general PDMPs lack 
a representation as a solution of a stochastic differential equation. Finally, Algorithm [2] obviously 
reduces to the SSA or a purely deterministic ODE method in the degenerate case if cither the 
PDMP is piecewise constant, i.e., / = const., or there are no jumps present at all, i.e., A = 0. 

The following theorem states conditions on the functions / and A, which are, together with a valid 
numerical method for solving (jlip . sufficient to guarantee almost sure convergence in the sense of 
Dcfmition[TJ The measure /i is assumed to be as specified in Section[2]and the properties discussed 
in Section [27X1 hold. 



Theorem 1. Let (A(i)) te [ 0i T] be a regular PDMP with phase space E. Let f be bounded, Lipschitz 
continuous and continuously differ entiable on E, i.e., there exist constants M,L such that for all 
(y,6), (z,ti) e E it holds 

\f(y,0)\ < M, (20) 
\f(y,6) -f(z,#) | < L\{y,6)-(z,$)\. (21) 

Further, assume that the jump rate A is bounded, bounded away from zero, Lipschitz continuous and 
continuously differ entiable, i.e., there exist constants A m i n , A max ,L such that for all (y,0), {z,$) £ 
E it holds 

< A m i n 

< X(y,0) < A max < oo, (22) 
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max{t„,t„} 



1 

max{t„ + i,i„ + i} 

■ 1 

tn+1 ~ in+1 I 



Figure 3: Even if the step size h is chosen small enough such that the post jump values of the 
piccewise constant component of the exact PDMP and its approximation coincide, we have that 
su Pte[ma X {t„J„},max{t„ +1 ,t-„ +1 }] W) ~ = Wn) - 0{t n+ i)\ > as in general \t n+1 -t n+1 \ > 
for all h > 0. 



\\(y,9)-\(z,#)\ < L\{y,9)-{z,d)\. (23) 

Then the algorithms [1] and [5] are well defined and the numerical approximation (-X"(i))te[o,T1 defined 
by Algorithm [2] converges almost surely to the exact PDMP constructed by Algorithm [T] for all 
initial values X(0) = xq G E in the sense of Definition [TJ Moreover, the continuous components 
of the PDMP and its approximation converge almost surely uniformly for h — > 0, i.e., 

lim sup \Y(t,u>) -Y(t,h,u)\ = (24) 
ft ~*°te[o,T] 

for almost all uj G fi. If, in addition, the continuous ODE method is of order p then also the 
asymptotic behaviour of ([T])-© and (|24p is of order p. 

We briefly comment on the nature of the conditions on / and A in the theorem. First of all, 
note that without loss of generality we can choose the same Lipschitz constant in (|2T|) and in 
(f23l) . Particularly this implies that the right hand side of (jTTJ) satisfies a Lipschitz condition 
with constant L, hence the global existence of a unique solution is guaranteed. Moreover, this 
solution is Lipschitz continuous with respect to the initial condition and, as /, A are continuously 
differcntiable, also diffcrcntiable with respect to the initial condition [12]. Secondly, the assumption 
that A is bounded away from zero is necessary in order for Step 2 in Algorithm [T] to be well-posed 
which presupposes, see |29| . on the one hand, that the derivative of the event function is non-zero 
at the event location, i.e., 

dr Jo 

On the other hand, it is further necessary for the right hand side of (fTTj) to be bounded in a 
ball around the event location in phase space, that is, in a ball around the point z G E x K + 
which the solution ip attains the moment its last component w hits the threshold. However, for 
the stochastic problem we are considering this event can occur at any point in phase space as 
the initial condition of the IVP ((TTj) as well as the threshold — \0gU2n-1 ar c random. Therefore 
the global boundedness of the right hand side of (fTT|) . i.e., the global boundedness of / and A, is 
required. 

We have already mentioned the connection to a convergence result for deterministic event detection 
|29| . Therein the authors prove that continuous methods conserve their order for problems where 
one jump event has to be approximated by event detection algorithms over the approximation 
interval. This result, however, is extended by Theorem Q] in two ways. First, random instead of 
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deterministic thresholds need to be considered for the simulation of PDMPs and secondly, jump 
times have to be calculated serially thus a simulation algorithm consists of a sequence of random 
event detection problems that take the result of the previous event detection as initial condition 
for the next one. Hence, for an error analysis we have to provide an analysis of the way the single 
errors of each event detection problem accumulate and are brought forward to a global level. 

4 Proof of Theorem [I] 

We begin by recalling some notation which is used throughout the proof. Let the finite interval 
[0, T] be fixed. We also fix an u> e f2 such that N = N(T,u>) < oo, i.e., the realised path 
{X(t, oj)) te [ .T] obtained by Algorithm [T] has only finitely many jumps in the interval [0,T]. By 
definition of a regular PDMP this includes all w except a potential set of measure zero. By 
t n (uj), n > 1, we denote the jump times of the exact PDMP. Further, for a fixed step size h 
of the continuous ODE method used to solve the system (fTTj) we denote by (X(t 7 h,oj)) te [o.T] 
the resulting path of the continuous approximation obtained by Algorithm [2] with jump times 
denoted by i ra (/i,w), n > 1. Finally, we denote by L^n-iXw), U2n{w), n > 1, the realisations of the 
i.i.d. standard uniform random variables defining the nth inter-jump time and the nth post-jump 
value, respectively. In the remainder of the proof we generally omit the dependency of any random 
variable on lu, as well as the dependency of the approximation on h. 

First, we briefly discuss the well-defincdncss of both Algorithms [1] and [5] are well-defined, i.e., 
we discuss the existence and uniqueness of trajectories for almost all uj G 51, which is inferred 
immediately from conditions ([2"T j) -(f23" |) . The Lipschitz conditions (|2"Tj) and (|23| guarantee the 
existence of a unique global solution ijj(t,x) to the IVP (fTTj) for all initial conditions x £ E and 
the existence of a numerical approximation ip(t,x) follows from the fact that the continuous ODE 
method is well-defined. Moreover, ([22]) implies that w(t, x) given by (fT0|) is strictly increasing with 
u>(0, x) = and w(t, x) — >• oo for t — > oo and hence (j4|) has a unique solution for any U2n-\ E (0, 1). 
For the corresponding property of the numerical approximation we note that a method may not 
produce continuous solutions that are increasing. However, the embedded discrete methods do 
produce strictly increasing and diverging iterates, which, e.g., immediately follows from the first 
order condition for Runge-Kutta methods [2], and condition (|22|) . By the continuity of w it follows 
that there exists a finite hitting time to every finite threshold. Uniqueness of this hitting time 
follows by the use of the generalised inverse in Step 2 of Algorithm^ Moreover, we remark that 
the same conditions imply that the approximate hitting times t n are each bounded in the step size 
h for a fixed realisation of the random threshold. 

We continue by explaining the basic idea of the proof before we present the arguments for the 
individual steps in detail. As usual for the convergence analysis of numerical methods we compare 
the approximation and the exact PDMP at discrete time points in [0,T], cf. the definition of the 
global error (|7|). Analogously to a standard approach in the convergence analysis we consider 
a 'local' error and show how the error is transported from the local to the global level. Then 
convergence follows from consistency and numerical stability of the method. However, the 'local' 
error we consider is not the error resulting over one discrete time step in the continuous ODE 
method but the error that arises over one interval between two successive jumps. This 'local' error 
corresponds to the global error of the ODE method over one inter-jump interval. 
Recall that we use a discretisation of the interval [0,T] with grid points given by the jump times 
t n and t n , n = 0, . . . , N, for the exact and the approximate PDMP, respectively. Therefore the 
time discretisation for the two processes we have to compare differ in the location of their grid 
points as well as possibly in the number of grid points, as in general t n ^ t n almost surely and also 
T < ijy may be possible with positive probability for some h > 0. The former is not a problem 
as we relate the nth inter-jump interval of the exact PDMP to the nth inter-jump interval of the 
approximation regardless of the interval endpoints, see Fig. [5] For this to be well defined the 
latter, however, is a problem as we thus need the same number of inter-jump intervals in each 
approximation as there are in the exact solution. Nevertheless, this problem is easily overcome, 
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as we can - without loss of generality - assume that Algorithm [5] is stopped only after it had at 
least N jumps. In this way we obtain an approximation with the same number of grid points and 
which is defined over an interval including [0,T] (at least for small enough h, cf. Section |4"~H) . 

The proof is organised as follows. In a first step of the proof we derive an error bound for the 
error in the approximation over one inter-jump interval, i.e., the 'local' error of our approximation 
method. To obtain such a bound we first need a bound on the error in the time until the next 
jump, i.e., a bound for 

|(£n+i — t n ) — {t n +i — t n ) \ . (25) 

Two distinct causes introduce errors in this approximation. On the one hand, the error due to the 
fact that we numerically approximate the hitting times t n +\ — t n by t n +i — t n and, on the other 
hand, due to perturbations in the initial conditions as in general X(t n ) ^ X(t n ). Therefore, we 
first consider in Section 14.11 consistency and numerical stability with respect to perturbations in 
the initial data for the hitting time approximation. Then, in Section 14.21 we consider the 'local' 
error of the method in the PDMP's phase space, whereas in Section POl we prove the first global 
convergence result, i.e., the limits Q and The convergence at the interval endpoint, i.e., flSJ), 
is proven in Section 14.41 Finally, in Section 14.51 we extend the discrete convergence result to a 
continuous result for the PDMP's continuous component, i.e., we prove the convergence (j24|) . 

4.1 Consistency of the hitting time approximation 

In this first part of the proof we consider (|25j) . i.e., the error in the time until the next jump of the 
exact PDMP started in x to the time until the next jump obtained by the numerical solution of (|11| 
started from a perturbed initial value x. Note that we consider the PDMP and its approximation 
just on one inter-jump interval and hence the components 9{t), 9(t) remain constant, however 
possibly distinct. 

To keep the presentation simple we omit the index n in the following, only to be reintroduced 
in Section 14.31 when we consider the global error of the method. Further, without loss of the 
generality we set for the exact and approximate inter-jump interval the left endpoint to and 
denote by t and t the jump time of the exact PDMP and its approximation, respectively. Hence 
the error ([25]) reduces to \t — t\ and the two times satisfy the equations 

— log £7 = w(t, x) and — logf/ = w(t, x) . 

That is, t and t are the respective hitting times of the random threshold — log U of the component 
w of the exact and approximate solution of the IVP (|lll) . 

To estimate \t — t\ we introduce the additional auxiliary time t as the time until the next jump of 
the exact process started in the perturbed initial value x which is given by 

— \ogU = w(t,x) . 

Then a simple application of the triangle inequality yields the initial estimate 

\t-t\ < \t-t\ + \t-t\. (26) 

Here the first term in the right hand side measures the error in the exact hitting time with respect 
to perturbations in the initial data whereas the second denotes the numerical error of the method 
without perturbations in the initial condition. We continue estimating the two terms separately 
and start with the latter. We introduce r(t, h, x) as the error of the continuous method's last 
component defined by 

w(t, x) = w(t, x) + r(t, h, x) 

and due to (fl~6| it is valid that \r(s, h,x)\ < err(t, h) for all s < t and all x € E. By definition of 
the hitting times and the Mean Value Theorem it follows that 

= w(t, x) + log U = w(t,x) -r(t,x, h) +log£7 

= w(t, x) + (t- t) Vtw(i9, x) - r(t, x, h) + log U 
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for some $ G [min{i, t}, max{t,t}]. As w(t,x) + \ogU = 0, this equality is equivalent to 

(t- t) V t w(ti, x) = -r(t, x, h), 
and we obtain due to the definition of w (TTU1) and condition (1221) the estimate 



,~ -> \r(t,x,h)\ errft.h) ,„„. 
\t - 1] < 1 v ;| < , ■ (27) 

Next we have to consider the term \t — t\ in the right hand side of (|2S|) . To bound this term we 
aim for an estimate on the distance in time by the distance of the initial conditions in phase space 
of pairs (i, x), (i,x) satisfying w(t, x) = w(t,x), i.e., an estimate of the form 

\t-t\ < Ci \x -x\ (28) 

for some positive random constant C\ depending only on t and U . Due to the assumptions 
in Theorem [1] the solution of the IVP (fTTj) depends differentiably on the initial value. Thus, 
proceeding analogously as for error estimate (|27|) we obtain by the Mean Value Theorem that 



= w(t, x) — w(t, x) 

= w(t,x) + V t w(^,C) (t- t) + V x w(tf,C) ■ (x-x) - w(t,x), 

for some mean values i? G [min{f, t}, max{<, t}], C G {y : y = sx + (1 — s)x, s G [0, 1]}. Then the 
Cauchy-Schwarz inequality yields 



t-t < 



(29) 



This gives an estimate of the form ([28)) if we bound the pre-factor in the right hand side of 
(|29l) uniformly. As the right hand side of the IVP (fTTj) is Lipschitz continuous in y and # with 
constant L it follows that x) depends Lipschitz continuously on the initial condition x with 
Lipschitz constant e Lt . Moreover, as by the assumptions of the theorem the solution is continuously 
diffcrcntiable with respect to the initial condition it follows that the derivatives satisfy 

\V x w(s,x)\ < c Ls Vie£, s>0. 

Hence we obtain for the pre-factor in the right hand side of (|2^|) the estimate 



V*u;0?,C) 



Kmc)) 



G 



L max{£,£} 



Am in 



Finally, it is easy to obtain an a priori bound for max{i, t}. Obviously, it is valid that max{i, t} < 
t + \t — t\ and as t, t are hitting times of the exact solution their difference is globally bounded 
due to (|22|) : A straightforward calculation shows that \t — t\ < A ^""' ~ A Amin (— log U) independently 
of the initial conditions x and x. These estimates applied to (f2T)| now yield 

~ e Lt 

\t-t\<- U~ s \x-x\, (30) 

A min 

where S = L A ^ lax ^ Amin and hence we have arrived at an estimate of the form (j28p . 

Overall an application of the estimates ([27|) and (f3Ul) to the right hand side of (|25|) yields the 
stability estimate 

\t-t\< — !— (c Lt U- s \x -x\+ err(t, h)) . (31) 
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4.2 The local error in phase space 

Next we derive a bound of the error over one inter-jump interval in the phase space E. We denote 
the post-jump values of the exact PDMP and its solution, started at initial conditions x and 
x, after one inter-jump interval by X(t) and X(t), respectively. Hence, given a realisation of a 
standard uniformly distributed random variable U' - here we use the notation U' to distinguish 
this random variable from U used defining the inter-jump time - the 'local' error in the phase 
space is given by 

\X(t) - = |0(t, x) + H{(f>{t, x), U') - M x) - H($(t, x), U')\ . (32) 

Note here H ( • , • ) is used as an abbreviation of the correct but cumbersome notation (0, H( • , ■)). 
To bound this error, we first estimate the difference in continuous components in (|32p . Thus, as 
a first estimate the triangle inequality yields 

\4>(t,x)-4>(t,x)\ < \(t>(t,x)-0,x)\ + \0,x)-$(t,x)\. (33) 

The second term in the right hand side of (|33|) is bounded by err(t, h) due to stability of the ODE 
method (flT)]) . Another application of the triangle inequality to the first term in the right hand 
side of (f3"3")) yields 

\4>{t,x) -</>(t,x) | < \<j>(t,x)- ( f>(t,x)\ + \ ( f>(t,x)-<l>(t,x)\. (34) 

As the solution of the IVP pip depends Lipschitz continuously on the initial condition the first 
term in the right hand side of (f3"4"|) satisfies 

\<f)(t,x) - 4>{t,x)\ < c Lt \x-x\. 

To estimate the second term in the right hand side of (|3"4"|) we employ the boundedness (|2T))) of / 
which yields 

\4>{t,x) - <j>(t,x)\ < M\t-t\. (35) 
Therefore we overall obtain for the left hand side of (f34|) the estimate 

\<t>{t,x) - 4>(t,x)\ < e Lt \x-x\+M\t-t\. 

Thus employing estimate (|3Tj) we get an error bound for the left hand side of (|33|) by 

m, x) - $(t, x)\ < e Lt (1 + X m ] n M U- s ) \x - x\ + (1 + X m ] n M) err(% h) , (36) 

which bounds the error due to the continuous part over one inter-jump interval. 

Finally, we consider the error introduced by the jump heights H for jumps starting at different 
points in phase space. In the following we do not bound the error due to the jump heights, but 
actually show that for sufficiently small step size this error vanishes. We denote for for all x, y S E 
the error by R(x, y, U') = \H(x,U') — H(y,U')\ which is finite almost surely and the considerations 
regarding the structure of the function H in Section 12.11 yield that 

R(<f>(t,x),<j>(t,x),U') = 

if \4>(t, x) — <f>(t, x) \ is smaller some bound depending on <fi(t, x) and U' . As we have already bound 
the difference \4>(t, x) — (j>(t, x)\, see 1|36|). we obtain that R(<j)(t, x), <p(t, x),V) = for small enough 
h and \x — x\. 

To conclude, overall we obtain for the 'local' error (f3"2"j) the estimate 

\X(t) - X(t)\ < c Lt (l + A m * n M U- s ) \x - x\ + (1 + \~l M) err(t, h) + R{cf>(t, x), $(% x), U') , 
where the last term vanishes for small enough h and \x — x\. 
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4.3 The discrete global error 

In order to prove convergence we have to be able to infer from the 'local' to the global error. Thus 
we now reintroduce the indices n, hence t n „ t n denote the nth jump times of the PDMP and its 
approximation, respectively, and X(t n ), X(t n ) denote the post-jump values of these processes at 
their jump times. We first show by induction over n that the global error at each jump time, 
i.e., \X(t n ) — X(t n )\ for all n = 1, . . . , N, converges to zero for h — > 0. Afterward we derive an 
upper bound for the global error uniform over all jump times to determine the asymptotic order 
of decrease of the global error. The splitting of the proof at this point into these two parts is 
necessary as in order to derive the asymptotic order of convergence we first need to show that the 
errors R due to different jump heights vanish for small enough step sizes h. However, this is a 
necessary condition of convergence and we can assume it holds for small enough step sizes once 
the convergence is established. 

For the induction basis we set X(0) = X(0), then (|3"T)) yields the estimate 

|X(ti)-^(?i)| < {l + \^M)err{t u h) + R{dp{t,x )M^ ),U 2 ). 

Obviously for h small enough the last term on the right hand side vanishes due to the convergence of 
the deterministic method. Then h — > implies that (l + A~?" n M) err(t\, h) — > as the approximate 
hitting time t\ = t\{h) is bounded in h. Therefore we obtain lim^o \X(ti) — X(ti)\ = . 

Next, for the induction step we assume that lim^o \X(t n ) — X(t n )\ = holds for some n > 1. 
Then, for small enough h we have that 

\X(t n+1 ) - X(t n+1 )\ < c^ 1 " 4 ") (1 + X m l M U^ +1 ) \X(t n ) - X(t n )\ + (1 + A m [ n M) err(t n+1 - t n , h) 
and hence lim,^ \X(t n+1 ) - X(t n+1 )\ = 0. 

Therefore overall we find that the method is convergent in the sense that 

lim max \X(t n ) - X(t n )\ = . (38) 

h— >0 n=l,...,N 

Analogously we show that also the exact and approximate jump times converge. The error estimate 
(jSTTl yields 



\tn+l — tn+l\ < — ^n) ~ {tn+1 ~ t n ) + t n — t n \ 

< e L(Wi-t») A -l n V -5 +i \ X (t n ) - X(t n )\ + \-\ n err(t n+1 - t n , h) + \t n - t n \. (39) 
Then p8p and another inductive argument yield that 

lim max \t n — t n \ = . (40) 

h->0n=l,...,N 



As these considerations leading to (|38|) and (|40|) arc valid for almost all u> S Vl this completes the 
proof of the properties (|7|) and ([9]). However, to satisfy the condition of pathwise convergence, 
Definition [l] it remains to show the convergence at the interval endpoint which is deferred 
to Section l4~4l At this point we continue to derive estimates for the global errors (0 and ([9]) to 
obtain a rate of convergence. 

To this end we recursively apply (|37[) to its right hand side arriving at 

\X(t n )-X(t n )\ < Z n \x -x \ (41) 

n 

+ Z n Zk 1 [(! + A mL M) err(t k - i fe _ 1; h) + R(X(t k _), X(t k -), U 2k ) 
fc=i 
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where for k = 1 . . . , N 

k 

Zu = c Ltfc H(l + \J n M U^_ x ) > 1 . (42) 
t=i 

As we are considering the asymptotic order of convergence for h — > we can assume that h is 
small enough such that R(X(t n -), X(t n -), t/2-n) = for all n = 1, . . . , N. Thus for xq = xq we 
obtain that for all 11 = 1, . . . , N 

n 

\X(t n )-X(t n )\ < (l + \- in M)Z n Y l Zk 1 err{t k -t k ^ 1 ,h), (43) 

k=l 

and therefore it holds that 

N 

max \X(t n ) -X(t n )\ < {l + \-] n M)Z N y / Z^ 1 err(t k -t k - 1 ,h). 

fc=l 

Then err(t, h) = 0{h,P) implies that 

max \X{t n )-X{t n )\ = 0{h p ). (44) 

n— 1,...,N 

Analogously we derive the order of convergence for the jump time approximations. Starting with 
to = t and recursively applying ([551) to its right hand side yields 

n 

\t n -t n \ < A-J n £[err(t fc -t k -x,h) + e^-^ U^_ x |A"(t fc _i) - X(t k . x )\ 

k=l 



Employing (|43[) to estimate the error in phase space yields for small enough h that for all n 
l,...,N 

n n—1 

\t n - t»\ < Yl KL [! + (! + A min M) ^ £ %i ^+1 e^"*')] err(t k t k . x ,h) . 

k=l j=k 

Denoting the pre- factor to err in each summand by W£ yields that for all n = 1, . . . , N 

N 



\tn -t n \<Y, W k err (tk ~ 4-1 , h) 



and thus err(t, h) = 0(h p ) implies 



fe=i 



max N \t n -t n \ = 0(h p ). (45) 



These considerations leading to (|4"4"| and (|4"5)) are valid for almost all ui G and thus (J7J) and 
follow with the proposed order. 



4.4 Convergence at the interval endpoint 

Finally for the method to satisfy the definition of pathwise convergence we establish convergence at 
the interval endpoint T. To this end we employ the estimate Q39p on the error for the approximation 
of the (N + l)th jump time. This yields 

\t N +i - t N +i\ < e L ( t2N + 1 - t ") A"^ U N S +1 \X(t N ) - X(t N )\ + X m l n err(t N+1 - t N , h) + \t N -t N \. 
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As all terms in the right hand side of this inequality converge to zero for h — > it follows that 
there exists an h* > 0, depending on Um+i, such that t^+i > T for all h < h*. In this case we 
then obtain for the error at the interval endpoint the estimate 

\X(T)-X(T)\ = \<f>(T-t N ,X(t N ))-$(T-t N ,X(ib N ))\ 

< \cf>(T-t N ,X(t N ))-ct>(T-t N ,X(t N ))\ + \<f>(T-t N ,X{t N ))-^{T-t N ,X(t N ))\ 

< M\t N -t N \+e L ( T - T ^\X(t N )-X(t N )\+err(T-t N ,h). (46) 

Here we have used in the last inequality on the first term an estimate analogous to (|35[) and on the 
second term the method's stability estimate (|16p . Due to previous results and the assumptions 
on the numerical method the final estimate (|4"6"]l converges to zero for h — > for almost all ui G O 
and if the method is of order p then it holds almost surely that 

\X(T) - X(T)\ = 0(h p ). 
4.5 The continuous global error 

We finish the proof extending the statement (|44[) to its uniform version (|24|) . Uniform convergence 
can only hold for the continuous component Y(t) and cannot hold for the discontinuous component 
6(t), cf. Fig. [31 which is also shown in a mathematically precise way below. 

To prove the results in this section we first derive a simple auxiliary result on the asymptotic 
behaviour of the approximate jump times. As stated in the last section the convergence of the 
jump times (|40|) implies that £jv > T for all h < h* for some h*. Moreover there exists an h** < h* 
such that for all h < h** it is valid that 

(i„-i,i„)n (t„_i,£0 ^0 Vn = l,...,iV. (47) 

That is, (|4"T)) says, that for small enough h the nth inter-jump interval of the exact PDMP and 
its approximation overlap for all n. We prove (|47[) by contradiction. Assume that there exists an 
n< N such that (i„_i,£„) fl (t n -i,t n ) = for all small h. Then it holds that either 

which implies that for all small h either 

< \t n - t n -i\ < \t n -i - t n -i\ or < |* n _i - t n \ < \t n - t„\. 

But, both are contradictions to lim; l _ ! .o niax n= o j ....Ar \t n — t n \ =0 and thus (|47[) holds. 

For the remainder of this section we assume that h < h**, hence ijv+i > T and (|47|) hold and all 
errors R vanish. Then, a partition of the time interval [0,T] is given by the discretisation points 

=: so < si := max{ii, t±} < . . . < := maxjiTv, tN} < sjv+i := T . 

It now immediately follows that a uniform convergence result cannot hold for the component 0(t): 
On every interval [s n , s n +i] it is valid that 

sup \9(t)-9(t)\ = \9(t n )-6(t n+1 )\>0 Vh<h**. 

ie[s„,s n+ i] 

as 0(s n ) = 6(s n ) and 9 and 9 jump exactly one time within the interval [s n ,s„+i] to the same 
post-jump value. By definition, one jumps at s n+ i and the other strictly before. Therefore their 
maximal difference over the interval [s n ,s n +i] is exactly the non-zero jump height independently 
ofh. 
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After these preliminary considerations we proceed to obtain the continuous convergence result (|24p . 
To this end we use the time grid so, ■ • ■ , sjv+i and estimate the uniform error of the continuous 
component Y(t) on each interval [s n ,s n+ i]. Note that on each of these intervals it holds that 
either 

i-J Ui+i < t n +i or ii.) t n+ i < t n+ i. 

We first consider case L). Recall that we use y(t,x), y(t,x) to denote the first component of the 
exact and numerical solution of the IVP (fTTj) and we use the symbol V to denote the maximum 
of two expressions. Then 

max \Y{t)-Y{t)\ = max \Y(t)-Y{t)\V max \Y(t)-Y(t)\ 

te[«„,a„+i] te[s n ,t n+1 ] tG[t„+i,t„+i] 

= max \y(s,X(s n )) - y(s,X(s n ))\ 
se[o,t„+i-s„] 

V max_ \y(s,X(t n+1 )) -y(s,X(t n+1 ))\. 

sG[0,i n+ i-t„+i] 

Both error terms in the right hand side can be estimated using the stability of the continuous 
ODE method (JT6]) which yields 

max \Y(t)-Y(t)\ < \e L ^-^ \X(s n ) - X(s n )\ 

Ve L (t „ +1 -?„ +l) \X(t n+1 )-X(t n+1 )\\ +err(s n+1 -s n ,h). 

As we assume that h is small enough such that 9(t n ) = 9(t n ) and hence 0(s n ) = 8(s n ) we obtain 
max \Y(t)-Y(t)\< L L ^^- S ^ \Y(s n ) - Y(s n )\ V 

i€[s„,s„ + i] 

e i(t„+i-fn+i)(| X (f n+1 ) -X(t n+1 )\ + \X(t n+1 ) - A(t n+1 )|)] +err(s n+1 -s n ,h) 
< e i ( s »+ 1 - s ")|F( S „)-Y( S „)| V 

e i( *»+ I_f »+ l) (M |t„+i -i„+i| + |*(i„+i) - X(t„+i)|)] + err(s„+i - s n , h). 
Due to the results ((44]) and (05j it follows that 

e £(t B+1 -f B+1 )( M |f n+i _ tn+i | + |jsf( tn+1 ) _ X(t n+1 )\) + err(s n+1 - s n , h) = 0(h p ) . 
Hence we obtain the recursive relation 

max \Y(t) - Y(t)\ < e L ^~ 3 ^ \Y(s n ) - Y(s n )\ + 0(h p ) . (48) 

te[s„,s„ + i] 

Secondly, for the case ii.), i.e., t n+ \ < t n+ \, we analogously obtain the estimate 
max \Y(t) - Y(t)\ < [ e L(s "+ 1 " s " ) \Y(s n ) - Y(s n )\ V 

i€[s„,s„ + i] 



e £(Wi-Wi) 



(\X(t n+1 ) - X(t n+1 )\ + \X(t n+1 ) - X(t n +i)\) + err(s n+1 - s„, h) 



To estimate the term \X{t n+ i) — X(t n+ i)\ we employ the Lipschitz condition p^|) . Hence we 
obtain again a recursive relation of the form (|48[) . Thus such a relation holds for the interval 
[s n ,s n +i] in either cases. 

Finally, recursively applying (|48|) to its right hand side yields the estimate 

max max \Y(t) - Y(t)\ < e LT \Y(0) - Y(0)\ + 0(h p ) . 

n=0,...,N te[s n ,s n +i] 

These considerations hold true for almost all u G hence the uniform convergence ([24")) follows 
for Y(0) = Y(0) and the proof of Theorem [T] is completed. 
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5 Extensions of the convergence theorem 



The results of Theorem Q] are easily generalised to arbitrary PDMPs with right continuous paths. 
That is, we now consider PDMPs that do not have any qualitative differences in their components: 
All components allow for discontinuities, i.e., fi is a Markov kernel onto R d+m - however still 
discretely supported -, and in between jumps the trajectories follow a deterministic motion given 
by an ODE 

'v\ = (h{y,of 
/ \Mv,oi 

i.e., the component 9 is not pieccwisc constant anymore. Clearly, a continuous convergence result 
such as cannot hold anymore. However, it is straightforward to see that the result in Theorem 
[T]on the asymptotic behaviour of the errors at the jump times, at the interval end point and the 
errors in the jump times is still valid if f\ and f 2 satisfy the conditions specified therein for /. 

In a further extension we consider the just discussed general right-continuous structure of a PDMP 
but assume that the distributions p,(x, •) of the jump heights are continuous on R d+m for all x £ E 
instead of discrete. Corollary 23.4 in [9] guarantees as in Section [27X1 the existence of a function 
H : E x (0, 1) — > M. d+m providing realisations of the random jump heights from realisations of 
standard uniformly distributed random variables. In order to derive a convergence result in this 
case we impose on H the following Lipschitz condition: there exists a deterministic constant 
Lh < oo such that for all x, y E M. d it holds almost surely that 

\H(x,U 2n ) - H(y,U 2n )\ < L H \x-y\ Vn>l. (49) 

Then by only a small change in the proof of Theorem [T] we are able to establish the following 
convergence result. 

Corollary 1. Let (A(i)) tg [ .T] and (X(t)) t€ [ x] be a regular PDMP and its approximation, re- 
spectively, which have right- continuous paths. The distribution of the jump heights is continuous 
and (|49[) is satisfied. Then, under the conditions of Theorem [JJ Algorithm [~J| converges pathwise 
in the sense of Definition^ i.e., the global error satisfies (O - ©. Moreover, if the embedded 
continuous ODE method is of order p, then the order of convergence is p. 



Proof. The proof of (JZJ)-© in case of this general class of PDMPs works analogously to the proof 
of Theorem [T] in Section The arguments employed become even less technical. In particular, 
the last paragraph in Section 14.21 dealing with the error in the jump heights becomes redundant 
as does the induction argument in Section 14.31 We restrict the presentation to point out the main 
difference. 

The proof proceeds as the proof in Section [4] with the main difference in estimating the local error 
in phase space, cf. Section 14.21 Using the Lipschitz condition (|49|) we first obtain from (|32j) the 
estimate 

\X(t)-X@)\ = (1 + L H ) \4>{t,x) - 4>(t,x)\ . 

Then, estimating the difference \(f>(t,x) — (f>(t,x)\ as in Section B~2l yields for the local error in phase 
space an estimate by 

\X(t) - X(t)\ < e Lt (l + L H )(1 + X~l M U- 5 ) \x - x\ + (1 + L H )(l + \ m { n M) err{% h) , 

cf. the local error estimate (|37|) . Starting with X(0) = X(0) and a recursive application of this 
last inequality, see Section 14. 3[ yields the global error bound 

N 

max \X{t n )-X{t n )\ < (l + L H ){l + X m l n M)Z N y Z' 1 err{t k -t k - U h). 

k=\ 

As this calculations are valid for almost all u) £ D, the convergence ([7]) follows and is of order p 
in case err(t, h) = 0(h p ). The proofs for the limits (|9]) and ([8]) work completely analogous as for 
Theorem HI see Sections 14.31 and 14741 respectively. □ 
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Figure 4: Coefficients and interpolation formulae of the continuous Runge-Kutta methods imple- 
mented for the numerical examples: (a) general Butcher tableau, (b) forward Eulcr method, (c) 
trapezoidal rule, (d) RadauIIa method, (e) LobattoIIIa method. 



6 Numerical examples 

To illustrate the theoretical findings on the order of convergence we have implemented simulation 
methods for PDMPs based on continuous Runge-Kutta methods of different order. For an IVP 

y = .f(t,y), y{0) = y 

with y(t), yo G M. d and t G [0, T] an s-stage continuous Runge-Kutta method is a discretisation 
scheme of the form 

8 

y((n+l)h) = y{nh) + h^2Pif(nh + Cih,ki), n = 0, . . . ,JV - 1, h = T/N (50) 

i=i 

with stage values fej, i = 1 . . . , s, given by 

s 

ki = y(nh) + h y ] a,j/(n/i + cjh, kj) (51) 
j'=i 

and combined with an interpolation formula 

s 

y(nh + £h) = y(nh) + h ^ 6^(0 f{nh + Cih, ha), < £_ < 1 (52) 
i=i 

for the approximation on the intervals between the discretisation points nh. The coefficients 
j3 = (j3i, . . . , f3 s ) T , A = (aij)i,j=i,...,a an d c = (ci,...,c s ) T are given by Butcher tableaus as 
in Fig. 2] (a). We implemented the explicit Eulcr method (order 1), the trapezoidal rule (order 
2), the 2-stage RadauIIa method (order 3) and the 3-stage LobattoIIIA method (order 4). The 
coefficients and interpolation polynomials for these methods are given in Fig.|3](b)-(e) taken 
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from [2:. Each experiment consists of the same trajectory simulated using Algorithm [2] based 
on the different continuous Rungc-Kutta methods with decreasing step sizes h. We compared 
the approximations to a reference solution as an exact solution is not available. The reference 
solution is an approximation of the same trajectory simulated with very high accuracy. We stopped 
decreasing the step size h for an implementation of Algorithm [2] when the approximations entered 
the error range of the reference solution. 

We have applied the methods to a stochastic hybrid version of the Hodgkin-Huxley model for 
a patch of neuronal membrane [H [24] . We consider the model for two different parameters sets 
denoted by (PI) and (P2). The numerical values of both parameter sets can be found in the 
appendix. The set (PI) is the original Hodgkin-Huxley model for the squid giant axon. Since its 
introduction in [15] this model is the standard example in neuroscience for models of excitable 
membranes. The second set (P2) is taken from |23j wherein the authors experimentally compare 
the performance of various pseudo-exact algorithms with respect to certain test statistics. 
The model describes the electrical potential difference across a neuronal membrane which is dynam- 
ically changed due to varying ionic currents through the membrane. These currents are mediated 
by membrane proteins that span across the membrane and in changing their shape these open 
and close pores. Such proteins are called ion channels and are named after the class of ions they 
allow selectively to pass. Changes in channel states occur in a random, memoryless fashion with 
instantaneous rates depending on the current potential difference. Conversely, the current through 
the membrane is proportional to the number of open states. For a more physiological explanation 
and derivation of the mathematical model we refer to [19] and [4] for its description using PDMPs. 
The hybrid version of the Hodgkin-Huxley neuron model is a 14-dimensional PDMP where a one- 
dimensional continuous variable Y(t) <E K models the difference in the electrical potential. The 
remaining, piecewise constant variables 9(t) £ R 13 record the states of the ion channels immersed 
in the membrane. In the Hodgkin-Huxley model there arc two different families of ion channels: 
sodium (Na + ) and potassium (K + ) channels with 8 and 5 distinct states, respectively. Their 
kinetic schemes are given in Fig. [6] in the appendix. Each component of 6(t) thus counts the 
number of channels in the specific state, i.e., 6i,...,6g correspond to the states of the Na + - 
channels and 9g, . . . , #13 to the states of the i4f + -channels. 

The characteristics of the PDMP are as follows. Firstly, the family of ODEs ([1]) defining the inter 
jump evolution of the PDMP is given by the equation 



Cy = -g Na Ps(y - % a ) - g K °i3(y - E K ) - .g L (y - e l) 

Secondly, the instantaneous jump rate A is given by 



I(t). 



(53) 
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Thirdly, we specify the transition measure fx. The point probability of the event that one channel 
changes from state i to state j is given by the transition rate of state i to state j times the number 
of channels in state i divided by the total instantaneous rate A. For example, the probability of 
the event of one channel changing from state 1 to 2 - conditional on the jump time being t - is 
given by 

3o m (r(f))fli(t-) rff Y Mo(t-)), {(-1, 1,0,..., of}) . 

All other events have zero probability, that is two channels do not change states simultaneously 
almost surely. 
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Figure 5: A first exemplary numerical experiment for a trajectory to the parameter set (PI). 
We have plotted the errors in phase space (top left), jump times (top right) and in the discrete 
component (bottom left). The path of the reference solution's continuous component calculated 
using MATLAB®'s ode45 is shown in the bottom right panel. 

We remark that (fSU)) differs from the general form ([1]) due to the added time dependent function 
I(t) which denotes the external current input to the system. However, for the numerical experi- 
ments we present the input is piecewise constant of the form I(t) = const. ■ I(Ti,t 2 ] (t) (monophasic 
input) and hence one may think of the model as PDMPs 'glued' together at the times T\, T 2 with 
the final state being the initial condition of the next. The reason for this initialising time is that 
in numerical experiments one wants to sample an initial condition from the membrane at rest, i.e., 
with I(t) = 0, for the start of the response to an input at time T\. 

We present two sample trajectories for the parameter set (PI). The reference solution for the tra- 
jectory in Fig. [5] was calculated using MATLAB®'s ode45 implementation with AbsTol / RelTol 
= 2.22045c -14 . We used the built-in event detection for the calculation of the hitting times. For 
the trajectory in Fig. [5] the reference solution was calculated using the LobattoIIIa method with 
step size h — 5c -6 . On the one hand the use of LobattoIIIA illustrates that high order methods 
yield good results in the absolute error range for reasonable equidistant step sizes. On the other 
hand the use of ode45, which employs automated step-size selection, illustrates that methods 
employing automated step size selection yield good accuracy results and thus can in principle be 
used. However, we found that for the same level of accuracy simulation times for the equidistant 
LobattolHA method were considerably shorter than for ode45. Thus an equidistant implementa- 
tion performs better as a simulation using MATLAB®'s time stepping algorithm contrary to the 
latter being specifically developed to speed up simulations. We note that these findings are not 
artefacts of the two specific paths presented but are persistent throughout all the simulations we 
conducted. 

The error plots in the upper panels of Fig. [SJ i.e., the error in phase space (J7J) (upper left panel) 
and the error in the jump times (|9]) (upper right panel) illustrate very clearly the theoretical result 
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Figure 6: A second exemplary numerical experiment for a trajectory to the parameter set (PI). 
We have plotted the errors in phase space (top left), jump times (top right) and in the discrete 
component (bottom left). The path of the reference solution's continuous component calculated 
using the presented PDMP method based on the continuous LobattoIIIA is shown in the bottom 
right panel. 

in Theorem Q] For step size h small enough such that the global error in the discrete component 
vanishes (lower left panel) the theoretical asymptotic order of convergence for the different methods 
is observed. For a guide of the eye we added lines of slope 1, 2, 3 and 4 beneath the errors of 
the methods with the respective orders. The plots of the second trajectory shown in Fig. [5] are 
a slightly less clear illustration of the asymptotic order of convergence. Here, the discrete error 
of the higher order methods has already vanished however reappears for a certain step size h. 
However, if we remove this outlier we observe the asymptotic order underlying the approximation 
as seen by comparison with the straight lines. 

As mentioned discussing the reference solution a naive use of automated step-size control with the 
intention of speeding up simulations and controlling the error can be misleading. However, it would 
be a task for an efficient and practicable pathwise step size detection and error control to detect 
such 'bad' step sizes as occur in the second trajectory and either avoid them or minimise their 
effect. These are crucial points for the implementation and performance of the algorithms. They 
demand for a further thorough investigation in this direction, e.g., an analysis of the shortcomings 
of the step-size detection as employed in ode45, which we have not attempted to do in this study. 
Finally, in a second example we consider a variant of the Hodgkin- Huxley model only dealing with 
currents due to A^a + -channels, that is the model reduces to 9 dimensions as 6g(t), . . . , Oxs(t) = 0. 
The channel density is higher and the current per channel smaller as for the first parameter set 
which overall renders the trajectories 'less noisy', cf. the sample trajectory of the reference solution 
in Figs. [5] and [6] with the one in Fig. [7] The reference solution was again calculated using the 
LobattoIIIA method with step size h = 5e~ 6 . Overall we find the same behaviour as in the first 
example. For small enough step sizes such that the errors in the discontinuous components vanish, 
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Figure 7: An exemplary numerical experiment for a trajectory to the parameter set (P2). We have 
plotted the errors in phase space (top left), jump times (top right) and in the discrete component 
(bottom left). The path of the reference solution's continuous component calculated using the 
presented PDMP method based on the continuous LobattoIIIA is shown in the bottom right 
panel. 



the predicted order of the asymptotic error behaviour is observed. However, note that even for the 
smallest step size considered the Euler method lacks the accuracy to approximate the trajectory 
correctly. Its error is orders of magnitude larger than the error of higher order Rungc-Kutta 
method with rather large step size. This observation strongly supports the use of higher order 
methods in practice. 
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Figure 8: Kinetic scheme of a (a) Aa + and (b) K + ion channel in the standard Hodgkin-Huxley 
model. The states 8 and 13 are the conducting states of channels, respectively. The rates a x = 
a x (y), b x = b x (y) are dependent on the transmembrane potential y. 

Appendix. 

The following are the rate functions and coefficients for the original Hodgkin-Huxley model of the 
squid giant axon taken from [15]. The rate functions are 

a n(y) = ioo(c( 10 -»)/ 10 -l) ' a m(y) = 10( - (25-b)/io_i) 7 a h(y) = 0.07c ^Z 20 , 

b n (y) = 0.125e-«/ 80 , b m (y) = 4c^/ 18 , b h (y) - 

and the constants used arc 



(PI) 



e (30- B )/10 +1 7 



£ , Na = H5mV, g Na = 4 pS, ?] Na = 300 ^m -2 , 

E K = -12 mV, g K = 18 pS, r] K = 30 /inT 2 , 

E L = mV, p L = 0.3 mS/cm 2 , 
C = 1 AiF/cm 2 . 



(PI) 



where r?Na and ryx denote the ./Va + and K + channel density, that is for the numerical example 
we considered a model with 300 iVa + -channels and 30 if + -channels. The input we consider is 
a monophasic current starting at t = 1 ms and lasting for 1 ms of strength 30 pA, i.e., I(t) = 
30-I (1)2] (t). 

The rate functions and parameters for the second example considered taken from [23j are given 
by 

/ \ _ 1.872(25.41-1/) i \ _ -0.549(27.74+1/) 

a m{y) — £,(25.41 — j/)/6.0G 2_ ' a h\V ) 1 _ c ( B + 27.74)/9.06 I 

, / x 3.973(21.001-;/) , / x 22,57 

u m\yj — 1 _ e ( v -21.001)/9.41 7 u h\yj — 1 + c (56-y)/12.5 ■ 

and the constants used arc 



(P2) 



^ Na = 144mV, # Na = 2.569 • 10~ 5 pS, ?7 Na = 1000 m™- 2 , 

£ , K =0mV, 7; K = 0pS, 7? K = 0/im" 2 , 

£ L = mV, 7? L = 1/(1953.49 • 10 3 ) mS/cm 2 , 1 ; 

C = 0.0714- 10" 6 /iF/cm 2 . 

Note that this examples considers Na + channels only. As input we consider a monophasic current 
starting at t = 0.1 ms and lasting for 0.1 ms of strength 35.1 • 10~ 6 pA, i.e., I(t) = 35.1 • 10~ 6 • 
1(0.1,0.2] (t). These are the parameters for the neuron considered in |23| . We note that the input 
current strength for the experiments is incorrectly reported in [23] and corrected in [3|. 
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